###########################################
# One week exposure 
###########################################

rm (list = ls())

library(Hmisc)
library(readstata13)
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(ggplot2)
library(did)
library(lfe)
library(gridExtra)

# load data
load("one_week_protests_riots_2024March1.RData")

##############
# Table A2
##############

mod_repression_indiv = felm(support_repression  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 + woman + age + education | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
tab_repression_indiv = mod_repression_indiv
mod_repression_indiv = tidy(mod_repression_indiv,conf.int = TRUE)
mod_repression_indiv

##############
# Table A3
##############

mod_workers_indiv = felm(support_workers  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 + woman + age + education | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
tab_workers_indiv = mod_workers_indiv
mod_workers_indiv = tidy(mod_workers_indiv,conf.int = TRUE)
mod_workers_indiv

##############
# Table A13
##############

mod_treat_predictors = felm(t_ind_after  ~ woman + age + education | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
mod_treat_predictors = tidy(mod_treat_predictors,conf.int = TRUE)
mod_treat_predictors

##############
# Table A14
##############

mod_eviction_indiv_nocontrols = felm(support_eviction  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
tab_eviction_indiv_nocontrols = mod_eviction_indiv_nocontrols
mod_eviction_indiv_nocontrols = tidy(mod_eviction_indiv_nocontrols,conf.int = TRUE)
mod_eviction_indiv_nocontrols

##############
# Table A15
##############

mod_students_indiv_nocontrols = felm(support_students  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
tab_students_indiv_nocontrols = mod_students_indiv_nocontrols
mod_students_indiv_nocontrols = tidy(mod_students_indiv_nocontrols,conf.int = TRUE)
mod_students_indiv_nocontrols

##############
# Table A16
##############

mod_eviction_comuna = felm(support_eviction  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 + woman + age + education | year_f + comuna | 0 | comuna, data = d, keepCX = TRUE)
tab_eviction_comuna = mod_eviction_comuna
mod_eviction_comuna = tidy(mod_eviction_comuna,conf.int = TRUE)
mod_eviction_comuna

##############
# Table A17
##############

mod_students_comuna = felm(support_students  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 + woman + age + education | year_f + comuna | 0 | comuna, data = d, keepCX = TRUE)
tab_students_comuna = mod_students_comuna
mod_students_comuna = tidy(mod_students_comuna,conf.int = TRUE)
mod_students_comuna

##############
# Table A18
##############

mod_eviction_indiv = felm(support_eviction  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 + woman + age + education | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
tab_eviction_indiv = mod_eviction_indiv
mod_eviction_indiv = tidy(mod_eviction_indiv,conf.int = TRUE)
mod_eviction_indiv

##############
# Table A19
##############

mod_students_indiv = felm(support_students  ~ year2016 + year2017 + year2019 + year2020 + year2021 + year2022 + woman + age + education | year_f + comuna | 0 | idencuesta, data = d, keepCX = TRUE)
tab_students_indiv = mod_students_indiv
mod_students_indiv = tidy(mod_students_indiv,conf.int = TRUE)
mod_students_indiv

##############
# Figure 1
##############

year = rep(2016:2022, each = 2)
t_ind = rep(0:1, 7)

# eviction
outcomes_eviction = mapply(function(y, t) mean(d$support_eviction[d$t_ind==t & d$year_f == y], na.rm = TRUE), year, t_ind)
outcomes_data_eviction = data.frame(outcomes_eviction,year,t_ind)
outcomes_data_eviction$groups = ifelse(outcomes_data_eviction$t_ind==0,"Control","Exposed")

fig_mean_eviction_by_year = ggplot(outcomes_data_eviction, aes(year,outcomes_eviction)) + 
  geom_line(aes(colour = groups), size = 1) + 
  geom_point(aes(colour = groups), size = 8) + 
  labs(x="", y="Average support (1-5 scale)") + 
  theme_bw() +  
  scale_fill_manual(values=c("black", "gray60") ) + 
  scale_colour_manual(name="Groups",  values =c("Control"="black", "Exposed"="gray60")) + 
  theme(legend.key = element_rect(colour = NA)) + 
  theme(text = element_text(size=40)) + 
  scale_x_continuous(breaks=c(2016:2022)) + 
  scale_y_continuous(limits = c(1, 3), breaks = c(1,1.5,2,2.5,3)) + 
  geom_vline(xintercept=2019, linetype="dashed", color = "black") + 
  xlab("Year") + 
  ggtitle("Support for violence from the state") + 
  theme(plot.title = element_text(hjust = 0.5))

# students
outcomes_students = mapply(function(y, t) mean(d$support_students[d$t_ind==t & d$year_f == y], na.rm = TRUE), year, t_ind)
outcomes_data_students = data.frame(outcomes_students,year,t_ind)
outcomes_data_students$groups = ifelse(outcomes_data_students$t_ind==0,"Control","Exposed")

fig_mean_students_by_year = ggplot(outcomes_data_students, aes(year,outcomes_students)) + 
  geom_line(aes(colour = groups), size = 1) + 
  geom_point(aes(colour = groups), size = 8) + 
  labs(x="", y="Average support (1-5 scale)") + 
  theme_bw() +  
  scale_fill_manual(values=c("black", "gray60") ) + 
  scale_colour_manual(name="Groups",  values =c("Control"="black", "Exposed"="gray60")) + 
  theme(legend.key = element_rect(colour = NA)) + 
  theme(text = element_text(size=40)) + 
  scale_x_continuous(breaks=c(2016:2022)) + 
  scale_y_continuous(limits = c(1, 3), breaks = c(1,1.5,2,2.5,3)) + 
  geom_vline(xintercept=2019, linetype="dashed", color = "black") + 
  xlab("Year") + 
  ggtitle("Support for violence from protesters") + 
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(fig_mean_eviction_by_year, fig_mean_students_by_year, nrow = 2)

##############
# Figure 2
##############

event_time = c(-3,-2,0,1,2,3)

# eviction
data_eviction = data.frame(event_time,
                           pe = mod_eviction_indiv$estimate[1:6],
                           se = mod_eviction_indiv$std.error[1:6],
                           ci_low = mod_eviction_indiv$conf.low[1:6],
                           ci_high = mod_eviction_indiv$conf.high[1:6])
data_eviction = rbind(data.frame(event_time = -1, pe=0, se=0, ci_low=0, ci_high=0), data_eviction)
data_eviction = data_eviction[order(data_eviction$event_time),]

fig_did_eviction <- ggplot(data = data_eviction,
                           mapping = aes(x = event_time, y = pe)) +
  geom_line(size = 0.5, alpha = 2, colour = "black") +
  theme_light()  +
  geom_vline(xintercept = 0, color = 'grey', size = 1.2, linetype = "dotted") + 
  geom_hline(yintercept = 0, colour="black",  linetype = "dotted")+
  geom_pointrange(aes(ymin = ci_low, ymax = ci_high), show.legend = FALSE, linetype= 1, size = 2,
                  color=ifelse(data_eviction$event_time<0,"gray60","black"))+
  theme(text = element_text(size=40)) + 
  xlab("Years to the social outbreak") +
  ylab("DiD Estimate") +
  scale_x_continuous(breaks = c(-3:4)) +
  scale_y_continuous(limits = c(-0.5,0.5)) + 
  ggtitle("Support for violence from the state") + 
  theme(plot.title = element_text(hjust = 0.5))

# students
data_students = data.frame(event_time,
                           pe = mod_students_indiv$estimate[1:6],
                           se = mod_students_indiv$std.error[1:6],
                           ci_low = mod_students_indiv$conf.low[1:6],
                           ci_high = mod_students_indiv$conf.high[1:6])
data_students = rbind(data.frame(event_time = -1, pe=0, se=0, ci_low=0, ci_high=0), data_students)
data_students = data_students[order(data_students$event_time),]

fig_did_students <- ggplot(data = data_students,
                           mapping = aes(x = event_time, y = pe)) +
  geom_line(size = 0.5, alpha = 2, colour = "black") +
  theme_light()  +
  geom_vline(xintercept = 0, color = 'grey', size = 1.2, linetype = "dotted") + 
  geom_hline(yintercept = 0, colour="black",  linetype = "dotted")+
  geom_pointrange(aes(ymin = ci_low, ymax = ci_high), show.legend = FALSE, linetype= 1, size = 2,
                  color=ifelse(data_students$event_time<0,"gray60","black"))+
  theme(text = element_text(size=40)) + 
  xlab("Years to the social outbreak") +
  ylab("DiD Estimate") +
  scale_x_continuous(breaks = c(-3:4)) +
  scale_y_continuous(limits = c(-0.5,0.5)) + 
  ggtitle("Support for violence from protesters") + 
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(fig_did_eviction, fig_did_students, nrow = 2)
